1 DATA

1.1 Pop sive over time dataset

data <- read.table(file=here::here("data", "BE_Census_Data.csv"), sep=",", header=TRUE)
  
head(data)
##    Block       Old_Label                  Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High1      High1              High
## 2 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High2      High2              High
## 3 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High3      High3              High
## 4 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High4      High4              High
## 5 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High5      High5              High
## 6 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High6      High6              High
##   Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1                  0               1     2/17/21         2/17/21       2/18/21
## 2                  0               1     2/17/21         2/17/21       2/18/21
## 3                  0               1     2/17/21         2/17/21       2/18/21
## 4                  0               1     2/17/21         2/17/21       2/18/21
## 5                  0               1     2/17/21         2/17/21       2/18/21
## 6                  0               1     2/17/21         2/17/21       2/18/21
##   Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1   DSC_0155       <NA>      NA      NA           101.9      NA      NA
## 2   DSC_0156       <NA>      NA      NA           104.7      NA      NA
## 3   DSC_0157       <NA>      NA      NA           105.8      NA      NA
## 4   DSC_0158       <NA>      NA      NA           101.9      NA      NA
## 5   DSC_0159       <NA>      NA      NA           102.2      NA      NA
## 6   DSC_0160       <NA>      NA      NA           101.7      NA      NA
##   HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1             100         NA          NA          no                     no
## 2             100         NA          NA          no                     no
## 3             100         NA          NA          no                     no
## 4             100         NA          NA          no                     no
## 5             100         NA          NA          no                     no
## 6             100         NA          NA          no                     no
##   Less_than_5
## 1          no
## 2          no
## 3          no
## 4          no
## 5          no
## 6          no
tail(data)
##      Block         Old_Label                  Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19      Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20      Med20
## 573 Block5   B5-Backup Fam A 8/25 BE B5 | G6 Med 21      Med21
## 574 Block5   B5-Backup Fam B 8/25 BE B5 | G6 Med 22      Med22
## 575 Block5   B5-Backup Fam C 8/25 BE B5 | G6 Med 23      Med23
## 576 Block5   B5-Backup Fam D 8/25 BE B5 | G6 Med 24      Med24
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571               Med                  5               6     8/25/21
## 572               Med                  5               6     8/25/21
## 573               Med                  5               6     8/25/21
## 574               Med                  5               6     8/25/21
## 575               Med                  5               6     8/25/21
## 576               Med                  5               6     8/25/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2  MC_Box1  MC_Box2
## 571         8/25/21       8/26/21   DSC_0273   DSC_0274  2.00000  2.00000
## 572         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 573         8/25/21       8/26/21   DSC_0277   DSC_0278 29.76978 32.92378
## 574         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 575         8/25/21       8/26/21   DSC_0267   DSC_0268 25.31021 46.10092
## 576         8/25/21       8/26/21   DSC_0283   DSC_0284 71.09684 55.37443
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571             4.0       2       2               4          8    0.500000
## 572             0.0      NA      NA              NA          0          NA
## 573            62.7      23      23              46         35    1.314286
## 574             0.0      NA      NA              NA         NA          NA
## 575            71.4      20      36              56         40    1.400000
## 576           126.5      76      57             133         41    3.243902
##     First_Throw Extinction_finalstatus Less_than_5
## 571          no                     no         yes
## 572     extinct                    yes         yes
## 573          no                     no          no
## 574     extinct                    yes         yes
## 575          no                     no          no
## 576          no                     no          no
dim(data)
## [1] 576  23
summary(data)
##     Block            Old_Label            Label            Population       
##  Length:576         Length:576         Length:576         Length:576        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##  Genetic_Diversity  Generation_Parents Generation_Eggs Date_Census       
##  Length:576         Min.   :0.0        Min.   :1.0     Length:576        
##  Class :character   1st Qu.:1.0        1st Qu.:2.0     Class :character  
##  Mode  :character   Median :2.5        Median :3.5     Mode  :character  
##                     Mean   :2.5        Mean   :3.5                       
##                     3rd Qu.:4.0        3rd Qu.:5.0                       
##                     Max.   :5.0        Max.   :6.0                       
##                                                                          
##  Date_Start_Eggs    Date_End_Eggs       Image_Box1         Image_Box2       
##  Length:576         Length:576         Length:576         Length:576        
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##     MC_Box1          MC_Box2       MC_Total_Adults     HC_Box1      
##  Min.   :  0.00   Min.   :  0.00   Min.   :  0.00   Min.   :  0.00  
##  1st Qu.: 26.52   1st Qu.: 25.99   1st Qu.: 51.05   1st Qu.: 22.50  
##  Median : 58.87   Median : 54.29   Median :101.50   Median : 57.00  
##  Mean   : 75.88   Mean   : 70.57   Mean   :128.37   Mean   : 72.23  
##  3rd Qu.:106.93   3rd Qu.: 99.40   3rd Qu.:170.25   3rd Qu.:101.00  
##  Max.   :327.40   Max.   :299.00   Max.   :514.40   Max.   :288.00  
##  NA's   :142      NA's   :141      NA's   :5        NA's   :145     
##     HC_Box2       HC_Total_Adults   Nb_parents     Growth_rate    
##  Min.   :  0.00   Min.   :  0.0   Min.   :  0.0   Min.   :0.0000  
##  1st Qu.: 25.00   1st Qu.: 63.0   1st Qu.: 68.0   1st Qu.:0.4379  
##  Median : 52.00   Median :100.0   Median :100.0   Median :0.9439  
##  Mean   : 68.08   Mean   :133.1   Mean   :138.6   Mean   :1.4639  
##  3rd Qu.: 96.75   3rd Qu.:169.0   3rd Qu.:189.0   3rd Qu.:2.3700  
##  Max.   :290.00   Max.   :499.0   Max.   :499.0   Max.   :6.8571  
##  NA's   :142      NA's   :48      NA's   :123     NA's   :147     
##  First_Throw        Extinction_finalstatus Less_than_5       
##  Length:576         Length:576             Length:576        
##  Class :character   Class :character       Class :character  
##  Mode  :character   Mode  :character       Mode  :character  
##                                                              
##                                                              
##                                                              
## 
#Remove populations killed due to the pathogen
data_kill<-data[data$Extinction_finalstatus=="kill",]
dim(data_kill)
## [1] 72 23
data<-data[data$Extinction_finalstatus!="kill",]

#Add pop growth rate
data$Nb_adults <- data$HC_Total_Adults
data$Lambda <- data$Nb_adults / data$Nb_parents
data$obs<-as.factor(1:nrow(data))
#summary(data)
#Remove 

#Check variables
str(data)
## 'data.frame':    504 obs. of  26 variables:
##  $ Block                 : chr  "Block3" "Block3" "Block3" "Block3" ...
##  $ Old_Label             : chr  "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
##  $ Label                 : chr  "BE B3 2/17 | G1  High1" "BE B3 2/17 | G1  High2" "BE B3 2/17 | G1  High3" "BE B3 2/17 | G1  High4" ...
##  $ Population            : chr  "High1" "High2" "High3" "High4" ...
##  $ Genetic_Diversity     : chr  "High" "High" "High" "High" ...
##  $ Generation_Parents    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Generation_Eggs       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Date_Census           : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_Start_Eggs       : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_End_Eggs         : chr  "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
##  $ Image_Box1            : chr  "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
##  $ Image_Box2            : chr  NA NA NA NA ...
##  $ MC_Box1               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Box2               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Total_Adults       : num  102 105 106 102 102 ...
##  $ HC_Box1               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Box2               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Total_Adults       : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Nb_parents            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Growth_rate           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ First_Throw           : chr  "no" "no" "no" "no" ...
##  $ Extinction_finalstatus: chr  "no" "no" "no" "no" ...
##  $ Less_than_5           : chr  "no" "no" "no" "no" ...
##  $ Nb_adults             : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Lambda                : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ obs                   : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
data$Genetic_Diversity <- as.factor(data$Genetic_Diversity)
data$Block <- as.factor(data$Block)
data$Population <- as.factor(data$Population)
data$Extinction_finalstatus <- as.factor(data$Extinction_finalstatus)


#Order levels 
data$Genetic_Diversity <- plyr::revalue(data$Genetic_Diversity, c("Med"="Medium"))
data$Genetic_Diversity <- factor(data$Genetic_Diversity, levels = c("High", "Medium", "Low"))
levels(data$Genetic_Diversity)
## [1] "High"   "Medium" "Low"
data<-droplevels(data) #remove previous levels 
str(data)
## 'data.frame':    504 obs. of  26 variables:
##  $ Block                 : Factor w/ 3 levels "Block3","Block4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Old_Label             : chr  "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
##  $ Label                 : chr  "BE B3 2/17 | G1  High1" "BE B3 2/17 | G1  High2" "BE B3 2/17 | G1  High3" "BE B3 2/17 | G1  High4" ...
##  $ Population            : Factor w/ 84 levels "High1","High13",..: 1 9 19 24 25 26 27 28 39 49 ...
##  $ Genetic_Diversity     : Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 3 3 3 ...
##  $ Generation_Parents    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Generation_Eggs       : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Date_Census           : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_Start_Eggs       : chr  "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
##  $ Date_End_Eggs         : chr  "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
##  $ Image_Box1            : chr  "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
##  $ Image_Box2            : chr  NA NA NA NA ...
##  $ MC_Box1               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Box2               : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ MC_Total_Adults       : num  102 105 106 102 102 ...
##  $ HC_Box1               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Box2               : int  NA NA NA NA NA NA NA NA NA NA ...
##  $ HC_Total_Adults       : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Nb_parents            : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ Growth_rate           : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ First_Throw           : chr  "no" "no" "no" "no" ...
##  $ Extinction_finalstatus: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Less_than_5           : chr  "no" "no" "no" "no" ...
##  $ Nb_adults             : int  100 100 100 100 100 100 100 100 100 100 ...
##  $ Lambda                : num  NA NA NA NA NA NA NA NA NA NA ...
##  $ obs                   : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
head(data)
##    Block       Old_Label                  Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High1      High1              High
## 2 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High2      High2              High
## 3 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High3      High3              High
## 4 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High4      High4              High
## 5 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High5      High5              High
## 6 Block3 B3_All_Red_Mix  BE B3 2/17 | G1  High6      High6              High
##   Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1                  0               1     2/17/21         2/17/21       2/18/21
## 2                  0               1     2/17/21         2/17/21       2/18/21
## 3                  0               1     2/17/21         2/17/21       2/18/21
## 4                  0               1     2/17/21         2/17/21       2/18/21
## 5                  0               1     2/17/21         2/17/21       2/18/21
## 6                  0               1     2/17/21         2/17/21       2/18/21
##   Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1   DSC_0155       <NA>      NA      NA           101.9      NA      NA
## 2   DSC_0156       <NA>      NA      NA           104.7      NA      NA
## 3   DSC_0157       <NA>      NA      NA           105.8      NA      NA
## 4   DSC_0158       <NA>      NA      NA           101.9      NA      NA
## 5   DSC_0159       <NA>      NA      NA           102.2      NA      NA
## 6   DSC_0160       <NA>      NA      NA           101.7      NA      NA
##   HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1             100         NA          NA          no                     no
## 2             100         NA          NA          no                     no
## 3             100         NA          NA          no                     no
## 4             100         NA          NA          no                     no
## 5             100         NA          NA          no                     no
## 6             100         NA          NA          no                     no
##   Less_than_5 Nb_adults Lambda obs
## 1          no       100     NA   1
## 2          no       100     NA   2
## 3          no       100     NA   3
## 4          no       100     NA   4
## 5          no       100     NA   5
## 6          no       100     NA   6
tail(data)
##      Block         Old_Label                  Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19      Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20      Med20
## 573 Block5   B5-Backup Fam A 8/25 BE B5 | G6 Med 21      Med21
## 574 Block5   B5-Backup Fam B 8/25 BE B5 | G6 Med 22      Med22
## 575 Block5   B5-Backup Fam C 8/25 BE B5 | G6 Med 23      Med23
## 576 Block5   B5-Backup Fam D 8/25 BE B5 | G6 Med 24      Med24
##     Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571            Medium                  5               6     8/25/21
## 572            Medium                  5               6     8/25/21
## 573            Medium                  5               6     8/25/21
## 574            Medium                  5               6     8/25/21
## 575            Medium                  5               6     8/25/21
## 576            Medium                  5               6     8/25/21
##     Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2  MC_Box1  MC_Box2
## 571         8/25/21       8/26/21   DSC_0273   DSC_0274  2.00000  2.00000
## 572         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 573         8/25/21       8/26/21   DSC_0277   DSC_0278 29.76978 32.92378
## 574         8/25/21       8/26/21       <NA>       <NA>       NA       NA
## 575         8/25/21       8/26/21   DSC_0267   DSC_0268 25.31021 46.10092
## 576         8/25/21       8/26/21   DSC_0283   DSC_0284 71.09684 55.37443
##     MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571             4.0       2       2               4          8    0.500000
## 572             0.0      NA      NA              NA          0          NA
## 573            62.7      23      23              46         35    1.314286
## 574             0.0      NA      NA              NA         NA          NA
## 575            71.4      20      36              56         40    1.400000
## 576           126.5      76      57             133         41    3.243902
##     First_Throw Extinction_finalstatus Less_than_5 Nb_adults   Lambda obs
## 571          no                     no         yes         4 0.500000 499
## 572     extinct                    yes         yes        NA       NA 500
## 573          no                     no          no        46 1.314286 501
## 574     extinct                    yes         yes        NA       NA 502
## 575          no                     no          no        56 1.400000 503
## 576          no                     no          no       133 3.243902 504

1.2 Remaining heterozygosity

#Upload growth rate end of experiment
data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
data_he$Population <- as.factor(data_he$Population)
data_he <- data_he[,c(1,3)]


### Additional: He remaining 
# New vector for the lost during exp
data_he$He_remain_exp <- NA

# He lost during exp
for(i in levels(data$Population)){
temp_data_pop <- subset(data, Population==i)
ifelse(sum(temp_data_pop$Generation_Parents)!=15,print("ERROR"),print(i))
temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_parents))
temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
data_he$He_remain_exp[data_he$Population==i] <- temp_he_remain_exp
rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
}
## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"
# Remaining He at the end of the experiment
data_he$He_remain_end <- data_he$He_remain * data_he$He_remain_exp

1.3 Phenotyping dataset

#Upload growth rate end of experiment
data_phenotyping <- read.table(file=here::here("data", "Adaptation_Data.csv"), sep=",", header=TRUE)


#Factors variables 
data_phenotyping$ID_Rep <- as.factor(data_phenotyping$ID_Rep)
data_phenotyping$Week <- as.factor(data_phenotyping$Week)
data_phenotyping$Block <- as.factor(data_phenotyping$Block)
data_phenotyping$Population <- as.factor(data_phenotyping$Population)
#Revalue
data_phenotyping$Genetic_Diversity <- as.factor(data_phenotyping$Divsersity)
data_phenotyping$Genetic_Diversity <-   plyr::revalue(data_phenotyping$Genetic_Diversity, c("Med"="Medium"))
data_phenotyping$Genetic_Diversity <- factor(data_phenotyping$Genetic_Diversity, levels = c("High", "Medium", "Low"))



#Add sum total 
data_phenotyping$Nb_ttx <- rowSums(data_phenotyping[,11:13], na.rm=TRUE)
#Proportion
data_phenotyping$p_adults <- data_phenotyping$Adults  / data_phenotyping$Nb_ttx
data_phenotyping$p_pupae <- data_phenotyping$Pupae / data_phenotyping$Nb_ttx
data_phenotyping$p_larvae <- data_phenotyping$Larvae / data_phenotyping$Nb_ttx

data_phenotyping_all <- data_phenotyping
data_phenotyping <- data_phenotyping[data_phenotyping$Week=="5-week",]

1.4 Merge dataset He

#Add heterozygosity
data <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)

#Add heterozygosity
data_phenotyping <- merge(x = data_phenotyping, y = data_he, by = "Population", all.x=TRUE)

1.5 Subset dataset generation

#Data beginning
data_G2 <- data[data$Generation_Eggs==2,]
#Data end 
data_G6 <- data[data$Generation_Eggs==6,]

1.6 Merge dataset phenotyping

## Merge dataset phenotyping
#merge dataset
temp_Start <- data_G2[,c("Block","Population",
                         "Genetic_Diversity","Lambda", "He_remain", "He_remain_end" )]
temp_Start$ID_Rep <- 1
temp_Start$Phenotyping <- "Initial"

temp_end <- data_phenotyping[,c("Block","Population",
                                "Genetic_Diversity","Lambda","ID_Rep", "He_remain", "He_remain_end" )]
temp_end$Phenotyping <- "Final"

data_evolution_Lambda <- rbind(temp_Start,temp_end)
data_evolution_Lambda$Phenotyping <- factor(data_evolution_Lambda$Phenotyping, 
                                            levels = c("Initial", "Final"))

2 PLOT

2.1 Population size

PLOT_Pop_size<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Population,
                                colour = Genetic_Diversity)) + 
  geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
  geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

# cowplot::save_plot(file = here::here("figures","OldPlot_PopulationSize.pdf"),PLOT_Pop_size, base_height = 10/cm(1),
#           base_width = 15/cm(1), dpi = 610)



SUM_Popsize <- Rmisc::summarySE(data, 
                        measurevar = c("Nb_adults"),
                       groupvars = c("Generation_Eggs",
                                     "Genetic_Diversity"), 
                       na.rm=TRUE)
SUM_Popsize
##    Generation_Eggs Genetic_Diversity  N Nb_adults        sd        se       ci
## 1                1              High 27 100.00000   0.00000  0.000000  0.00000
## 2                1            Medium 23 100.00000   0.00000  0.000000  0.00000
## 3                1               Low 34 100.00000   0.00000  0.000000  0.00000
## 4                2              High 27 344.29630  73.77294 14.197610 29.18360
## 5                2            Medium 23 282.04348 100.75735 21.009360 43.57075
## 6                2               Low 34 282.61765  78.53006 13.467794 27.40043
## 7                3              High 27 151.85185  80.96899 15.582489 32.03027
## 8                3            Medium 23 118.73913  88.54491 18.462891 38.28969
## 9                3               Low 34 104.61765  85.13341 14.600260 29.70445
## 10               4              High 26 114.00000  80.20224 15.728954 32.39439
## 11               4            Medium 23  62.65217  64.50483 13.450188 27.89398
## 12               4               Low 34  46.38235  39.63241  6.796903 13.82840
## 13               5              High 27 103.62963  51.56486  9.923661 20.39838
## 14               5            Medium 19  64.68421  49.88437 11.444259 24.04350
## 15               5               Low 30  53.16667  45.98207  8.395139 17.16999
## 16               6              High 27 147.66667  42.60282  8.198916 16.85311
## 17               6            Medium 18  73.61111  44.05496 10.383855 21.90802
## 18               6               Low 27  79.88889  44.43520  8.551559 17.57798
PLOT_Pop_size_average <- PLOT_Pop_size + 
  geom_point(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 5, position = position_dodge(0.4)) +
  geom_errorbar(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                ymin = Nb_adults-ci, ymax = Nb_adults+ci, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 1, width=.2,  position = position_dodge(0.4)) +
  geom_line(data = SUM_Popsize, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
           linetype = "dashed", size = 1.5, position = position_dodge(0.4)) 
PLOT_Pop_size_average

# cowplot::save_plot(file = here::here("figures","OldPlot_PopulationSize_average.pdf"),   PLOT_Pop_size_average, base_height = 10/cm(1),           base_width = 15/cm(1), dpi = 610)

2.2 Extinction

## Extinction yes no
PLOT_Extinction<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs), 
                                y = Nb_adults, 
                                group = Population,
                                colour = Extinction_finalstatus)) + 
  facet_wrap(~Genetic_Diversity, ncol=3) + 
  #geom_point(position = position_dodge(0.4), size = 0.4, alpha = 0.7) +
  geom_line(position = position_dodge(0.1), size = 0.4, alpha = 0.8) +
  ylab("Population size") +
  scale_color_manual(values = c("black", "red")) +
  xlab("Generation") +
  theme(legend.position = "none") +
  theme_LO
PLOT_Extinction

# #
# cowplot::save_plot(here::here("figures","OldPlot_Extinction.pdf"),   PLOT_Extinction, base_height = 6/cm(1),
#           base_width = 20/cm(1), dpi = 610)
# 
# #

2.3 Growth rate

tapply(data_G2$Lambda,data_G2$Genetic_Diversity,mean)
##     High   Medium      Low 
## 3.442963 2.820435 2.826176
PLOT_Growth<-ggplot2::ggplot(data_G2, aes(x = Genetic_Diversity, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  geom_boxplot(aes(group = Genetic_Diversity),outlier.colour = NA) +
  geom_jitter(width = 0.25, size = 1, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Genetic diversity") +
  theme_LO
PLOT_Growth

# 
# cowplot::save_plot(here::here("figures","OldPlot_Growth_Start.pdf"),   PLOT_Growth, base_height = 10/cm(1),
#           base_width = 8/cm(1), dpi = 610)

2.4 Growth rate evolution

SUM_MEAN_start_end <- Rmisc::summarySE(data_evolution_Lambda, measurevar = c("Lambda"),
                       groupvars = c("Genetic_Diversity","Phenotyping"), 
                       na.rm=TRUE)
SUM_POP_start_end <- Rmisc::summarySE(data_evolution_Lambda, measurevar = c("Lambda"),
                       groupvars = c("Genetic_Diversity","Phenotyping",
                                     "Population"), 
                       na.rm=TRUE)





PLOT_Growth_startend<-ggplot(SUM_POP_start_end, aes(x = Phenotyping, 
                                 y = Lambda, 
                                 group = Population,
                                 colour = Genetic_Diversity)) + 
  geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
  geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
  ylab("Growth rate") +
  xlab("Phenotyping") +
  theme_LO
PLOT_Growth_startend

PLOT_Growth_startend_average <- PLOT_Growth_startend + 
  geom_point(data = SUM_MEAN_start_end, aes(x = Phenotyping, 
                                y = Lambda, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 5, position = position_dodge(0.4)) +
  geom_errorbar(data = SUM_MEAN_start_end, aes(x = Phenotyping, 
                                ymin = Lambda-ci, ymax = Lambda+ci, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
             size = 1, width=.2,  position = position_dodge(0.4)) +
  geom_line(data = SUM_MEAN_start_end, aes(x = Phenotyping, 
                                y = Lambda, 
                                group = Genetic_Diversity,
                                colour = Genetic_Diversity), 
           linetype = "dashed", size = 1.5, position = position_dodge(0.4)) 
PLOT_Growth_startend_average

# 
# cowplot::save_plot(here::here("figures","OldPlot_Lambda_Evolution.pdf"),
#           PLOT_Growth_startend_average, base_height = 10/cm(1),
#           base_width = 15/cm(1), dpi = 610)

2.5 Life stage proportion

SUM_Prop_adults<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_adults"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_adults
##   Genetic_Diversity   Week   N  p_adults         sd          se         ci
## 1              High 4-week  50 0.6508098 0.17776618 0.025139934 0.05052059
## 2              High 5-week 105 0.9632040 0.05259208 0.005132461 0.01017786
## 3            Medium 4-week  24 0.7121270 0.15722597 0.032093616 0.06639070
## 4            Medium 5-week  49 0.9563136 0.06160975 0.008801393 0.01769639
## 5               Low 4-week  30 0.6568606 0.19029917 0.034743716 0.07105888
## 6               Low 5-week  59 0.9582021 0.07845154 0.010213520 0.02044458
SUM_Prop_pupae<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_pupae"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_pupae
##   Genetic_Diversity   Week   N    p_pupae         sd          se          ci
## 1              High 4-week  50 0.28825555 0.15629878 0.022103985 0.044419621
## 2              High 5-week 105 0.01458769 0.02407445 0.002349426 0.004658999
## 3            Medium 4-week  24 0.23019148 0.14506173 0.029610601 0.061254195
## 4            Medium 5-week  49 0.02283167 0.03485130 0.004978757 0.010010462
## 5               Low 4-week  30 0.29380227 0.17562943 0.032065401 0.065581108
## 6               Low 5-week  59 0.01675308 0.02459640 0.003202178 0.006409856
SUM_Prop_larvae<- Rmisc::summarySE(data_phenotyping_all, 
                        measurevar = c("p_larvae"),
                       groupvars = c("Genetic_Diversity","Week"), 
                       na.rm = TRUE)
SUM_Prop_larvae
##   Genetic_Diversity   Week   N   p_larvae         sd          se          ci
## 1              High 4-week  50 0.06093466 0.04274850 0.006045551 0.012148989
## 2              High 5-week 105 0.02220828 0.04324043 0.004219833 0.008368088
## 3            Medium 4-week  24 0.05768150 0.05999020 0.012245449 0.025331640
## 4            Medium 5-week  49 0.02085474 0.04462923 0.006375605 0.012819012
## 5               Low 4-week  30 0.04933713 0.08198954 0.014969174 0.030615398
## 6               Low 5-week  59 0.02504479 0.06691319 0.008711355 0.017437671
SUM_Prop <- data.frame(SUM_Prop_adults[,1:4],
                       p_pupae=SUM_Prop_pupae$p_pupae,
                       p_larvae=SUM_Prop_larvae$p_larvae)


SUM_Prop<-tidyr::gather(SUM_Prop,"Stage", "Proportion",4:6,-"Genetic_Diversity")
SUM_Prop
##    Genetic_Diversity   Week   N    Stage Proportion
## 1               High 4-week  50 p_adults 0.65080979
## 2               High 5-week 105 p_adults 0.96320403
## 3             Medium 4-week  24 p_adults 0.71212702
## 4             Medium 5-week  49 p_adults 0.95631359
## 5                Low 4-week  30 p_adults 0.65686059
## 6                Low 5-week  59 p_adults 0.95820213
## 7               High 4-week  50  p_pupae 0.28825555
## 8               High 5-week 105  p_pupae 0.01458769
## 9             Medium 4-week  24  p_pupae 0.23019148
## 10            Medium 5-week  49  p_pupae 0.02283167
## 11               Low 4-week  30  p_pupae 0.29380227
## 12               Low 5-week  59  p_pupae 0.01675308
## 13              High 4-week  50 p_larvae 0.06093466
## 14              High 5-week 105 p_larvae 0.02220828
## 15            Medium 4-week  24 p_larvae 0.05768150
## 16            Medium 5-week  49 p_larvae 0.02085474
## 17               Low 4-week  30 p_larvae 0.04933713
## 18               Low 5-week  59 p_larvae 0.02504479
SUM_Prop$Stage <- factor(SUM_Prop$Stage, levels = c("p_adults", "p_pupae", "p_larvae"))


PLOT_prop <- ggplot2::ggplot(data = SUM_Prop, aes(x = Genetic_Diversity, 
                                                  y = Proportion, fill = Stage)) +
  facet_wrap(~ Week, ncol=2) +
  geom_bar(stat="identity") + 
  xlab("Level of genetic diversity") +
  labs(color="Stage of individuals") +
  ylab("Proportion of each stage") +
  scale_fill_manual(values= c("#4EABB9", "#E95C2B", "#CBA938"), 
                    breaks = c("p_adults", "p_pupae", "p_larvae"),
                    labels = c( "Adult", "Pupae","Larvae")) + 
  theme_LO
PLOT_prop

# cowplot::save_plot(here::here("figures","OldPlot_Life_stage_proportion.pdf"),   PLOT_prop,
#           base_height = 8/cm(1), base_width = 16/cm(1), dpi = 610)
# 

2.6 Remaining heterozygosity

# All but with: He beginning for both and with pseudoreplication for final
ggplot2::ggplot(data_evolution_Lambda, aes(x = He_remain, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  facet_wrap(~ Phenotyping, ncol = 1)+
  geom_point(size = 2, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Remaining heterozygosity") +
  theme_LO

PLOT_He_Initial <- ggplot2::ggplot(data_G2, aes(x = He_remain, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  geom_point(size = 2, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Remaining heterozygosity") +
  theme_LO
PLOT_He_Initial

# Remove pseudoreplication: average of each population
SUM_POP_He_Mean<- Rmisc::summarySE(data_evolution_Lambda[data_evolution_Lambda$Phenotyping=="Final",],
                            measurevar = c("Lambda"),
                       groupvars = c("Genetic_Diversity","He_remain_end",
                                     "Population"), 
                       na.rm=TRUE)

PLOT_He_Final <- ggplot2::ggplot(SUM_POP_He_Mean, aes(x = He_remain_end, y = Lambda, 
                         colour = Genetic_Diversity)) + 
  geom_point(size = 2, alpha = 0.8) +
  ylab("Growth rate") +
  theme(legend.position = "none") +
  xlab("Remaining heterozygosity") +
  theme_LO
PLOT_He_Final

# ALL WITHOUT PSEUDO
SUM_POP_He_ALL<- Rmisc::summarySE(data_evolution_Lambda,
                            measurevar = c("Lambda"),
                       groupvars = c("Genetic_Diversity","Phenotyping", "He_remain", "He_remain_end",
                                     "Population"), 
                       na.rm=TRUE)

SUM_POP_He_ALL$HE <- ifelse(SUM_POP_He_ALL$Phenotyping=="Initial",SUM_POP_He_ALL$He_remain,SUM_POP_He_ALL$He_remain_end)


PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda, 
                         colour = Genetic_Diversity, 
                         shape = Phenotyping)) + 
  geom_point(size = 3, alpha = 0.8) +
  scale_shape_manual(values = c(1, 19)) + 
  ylab("Growth rate") +
  xlab("Remaining heterozygosity") +
  theme_LO
PLOT_He_ALL

# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",], aes(x = HE, y = Lambda,
#                          colour = Genetic_Diversity,
#                          shape = Phenotyping)) +
#   geom_point(size = 3, alpha = 0.8) +
#   scale_shape_manual(values = c(19)) +
#   ylab("Growth rate") +
#   xlab("Remaining heterozygosity") +
#   theme_LO
# PLOT_He_ALL

3 ANALYSIS

3.1 Probability of extinction

############ 
###### Clean dataset
############ 
#Prepare clean dataset
data_proba_extinction <- data_G6[,c("Block","Population","Genetic_Diversity","Extinction_finalstatus")]
data_proba_extinction$Extinction_finalstatus <- as.factor(data_proba_extinction$Extinction_finalstatus)
data_proba_extinction$y <- ifelse(data_proba_extinction$Extinction_finalstatus == "yes", 1, 0)



############ 
###### Analysis
############ 
#Analysis
mod0 <- lme4::glmer(y ~ Genetic_Diversity + Block + (1|Population),  
      data = data_proba_extinction, family = binomial)


summary(mod0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: y ~ Genetic_Diversity + Block + (1 | Population)
##    Data: data_proba_extinction
## 
##      AIC      BIC   logLik deviance df.resid 
##     58.8     73.4    -23.4     46.8       78 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1010 -0.3208 -0.2215  0.0000  4.5137 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Population (Intercept) 0        0       
## Number of obs: 84, groups:  Population, 84
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)   
## (Intercept)              -22.3274  5077.2816  -0.004  0.99649   
## Genetic_DiversityMedium   19.3132  5077.2815   0.004  0.99696   
## Genetic_DiversityLow      19.5192  5077.2815   0.004  0.99693   
## BlockBlock4                0.7402     1.2714   0.582  0.56046   
## BlockBlock5                3.0006     1.1265   2.664  0.00773 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM -1.000                     
## Gntc_DvrstL -1.000  1.000              
## BlockBlock4  0.000  0.000  0.000       
## BlockBlock5  0.000  0.000  0.000  0.737
## convergence code: 0
## boundary (singular) fit: see ?isSingular
#Test genetic diversity effect
mod1 <- lme4::glmer(y ~ Block + (1|Population),  
      data = data_proba_extinction, family = binomial)
anova(mod0, mod1, test = "Chisq")
## Data: data_proba_extinction
## Models:
## mod1: y ~ Block + (1 | Population)
## mod0: y ~ Genetic_Diversity + Block + (1 | Population)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## mod1    4 66.903 76.626 -29.451   58.903                        
## mod0    6 58.833 73.418 -23.417   46.833 12.069  2   0.002394 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#We keep genetic diversity 


mod0 <- lme4::glmer(y ~ Genetic_Diversity-1 + Block + (1|Population),  
      data = data_proba_extinction, family = binomial)
summary(mod0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: y ~ Genetic_Diversity - 1 + Block + (1 | Population)
##    Data: data_proba_extinction
## 
##      AIC      BIC   logLik deviance df.resid 
##     58.8     73.4    -23.4     46.8       78 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.1010 -0.3208 -0.2215  0.0000  4.5137 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Population (Intercept) 5.55e-15 7.45e-08
## Number of obs: 84, groups:  Population, 84
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)   
## Genetic_DiversityHigh   -3.737e+01  8.535e+06   0.000  1.00000   
## Genetic_DiversityMedium -3.014e+00  1.129e+00  -2.669  0.00760 **
## Genetic_DiversityLow    -2.808e+00  1.066e+00  -2.635  0.00843 **
## BlockBlock4              7.402e-01  1.271e+00   0.582  0.56046   
## BlockBlock5              3.001e+00  1.127e+00   2.664  0.00773 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             Gnt_DH Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM  0.000                     
## Gntc_DvrstL  0.000  0.767              
## BlockBlock4  0.000 -0.724 -0.786       
## BlockBlock5  0.000 -0.843 -0.871  0.737
## convergence code: 0
## boundary (singular) fit: see ?isSingular
############ 
###### Predicted value
############ 
#Extract probability of extinction 
(p_high <- plogis(-3.950e+01))
## [1] 7.004352e-18
(p_med <- plogis(-3.014e+00))
## [1] 0.04679739
(p_low <- plogis(-2.808e+00))
## [1] 0.0568934
predict(mod0, type = "response", 
        re.form = NA,
        newdata = data.frame(Genetic_Diversity = levels(data_proba_extinction$Genetic_Diversity), 
                             Block = "Block4"))
##            1            2            3 
## 2.220446e-16 9.329490e-02 1.122446e-01
############ 
###### Posthoc
############ 
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean      SE  df asymp.LCL asymp.UCL
##  High              -36.12 8534537 Inf -20378273  20378200
##  Medium             -1.77       1 Inf        -3         0
##  Low                -1.56       1 Inf        -3         0
## 
## Results are averaged over the levels of: Block 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate      SE  df z.ratio p.value
##  High - Medium  -34.352 8534537 Inf  0.000  1.0000 
##  High - Low     -34.558 8534537 Inf  0.000  1.0000 
##  Medium - Low    -0.206       1 Inf -0.274  0.9594 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log odds ratio (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.2 Population size G6

############ 
###### Analysis
############ 
#Analysis
# mod0 <- lme4::glmer(Nb_adults ~ Genetic_Diversity + Block + (1|Population),  
#       data = data_G6, family = "poisson")


mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block, 
             data = data_G6[data_G6$Extinction_finalstatus == "no",])


mod1 <- lm(log(Nb_adults) ~ Block , 
             data = data_G6[data_G6$Extinction_finalstatus == "no",])

anova(mod0, mod1) # Effect of genetic diversity 
## Analysis of Variance Table
## 
## Model 1: log(Nb_adults) ~ Genetic_Diversity + Block
## Model 2: log(Nb_adults) ~ Block
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     66 29.706                                  
## 2     68 42.118 -2   -12.412 13.789 9.914e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block, 
             data = data_G6[data_G6$Extinction_finalstatus == "no",])


summary(mod0)
## 
## Call:
## lm(formula = log(Nb_adults) ~ Genetic_Diversity + Block, data = data_G6[data_G6$Extinction_finalstatus == 
##     "no", ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.2417 -0.2947  0.1088  0.4381  1.2624 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              5.11532    0.17980  28.451  < 2e-16 ***
## Genetic_DiversityMedium -0.94813    0.20540  -4.616 1.86e-05 ***
## Genetic_DiversityLow    -0.80532    0.18719  -4.302 5.72e-05 ***
## BlockBlock4             -0.02077    0.18447  -0.113   0.9107    
## BlockBlock5             -0.53922    0.21414  -2.518   0.0142 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6709 on 66 degrees of freedom
## Multiple R-squared:  0.3362, Adjusted R-squared:  0.2959 
## F-statistic: 8.356 on 4 and 66 DF,  p-value: 1.623e-05
############ 
###### Posthoc
############ 
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE df lower.CL upper.CL
##  High                4.93 0.130 66     4.61     5.25
##  Medium              3.98 0.159 66     3.59     4.37
##  Low                 4.12 0.136 66     3.79     4.46
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE df t.ratio p.value
##  High - Medium    0.948 0.205 66  4.616  0.0001 
##  High - Low       0.805 0.187 66  4.302  0.0002 
##  Medium - Low    -0.143 0.207 66 -0.690  0.7704 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.3 Population size over time

PLOT_Pop_size_average

fit <- glm(Nb_adults ~ Generation_Eggs, data = data, family = "poisson")
#second degree
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), data = data, family = "poisson")
#third degree
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE), data = data, family = "poisson")
#exp
fit6 <- glm(Nb_adults~exp(Generation_Eggs), data = data, family = "poisson")
#log
fit7 <- glm(Nb_adults~log(Generation_Eggs), data = data, family = "poisson")


#generate range of 50 numbers starting from 30 and ending at 160
xx <- seq(1,6, length=100)
plot(data$Generation_Eggs,
     data$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")

# Check the r squared
summary(fit)$adj.r.squared
## NULL
summary(fit2)$adj.r.squared
## NULL
summary(fit3)$adj.r.squared
## NULL
summary(fit4)$adj.r.squared # best r square
## NULL
summary(fit5)$adj.r.squared 
## NULL
summary(fit6)$adj.r.squared 
## NULL
summary(fit7)$adj.r.squared 
## NULL
#
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")



#If we dont consider Gen=1 
PLOT_Pop_size_average

fit <- lm(log(Nb_adults+1) ~ Generation_Eggs, data = data[data$Generation_Eggs>=2,])
#second degree
fit2 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#third degree
fit3 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,3,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit4 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,4,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit5 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,5,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#exp
fit6 <- lm(log(Nb_adults+1)~exp(Generation_Eggs), data = data[data$Generation_Eggs>=2,])
#log
fit7 <- lm(log(Nb_adults+1)~log(Generation_Eggs), data = data[data$Generation_Eggs>=2,])
#generate range of 50 numbers starting from 30 and ending at 160
xx <- seq(1,6, length=100)
plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
     log(data[data$Generation_Eggs>=2,]$Nb_adults+1),pch=19,ylim=c(0,8))
lines(xx, predict(fit, data.frame(Generation_Eggs=xx)), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx)), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx)), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx)), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx)), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx)), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx)), col="pink")

# Check the r squared
summary(fit)$adj.r.squared
## [1] 0.1281606
summary(fit2)$adj.r.squared
## [1] 0.3033469
summary(fit3)$adj.r.squared
## [1] 0.3017434
summary(fit4)$adj.r.squared # best r square
## [1] 0.3031656
summary(fit5)$adj.r.squared 
## [1] 0.3031656
summary(fit6)$adj.r.squared 
## [1] 0.01432146
summary(fit7)$adj.r.squared 
## [1] 0.1763784
# data$gen_square <- data$Generation_Eggs^2
# fit2 <- lm(log(Nb_adults+1)~poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs + gen_square, data = data[data$Generation_Eggs>=2,])
# #Add the interaction with Genetic diversity 
# fit2 <- glm(Nb_adults~ Generation_Eggs*Genetic_Diversity + 
#              gen_square*Genetic_Diversity, 
#            data = data, family = "poisson")
# 
# 
# summary(fit2)


####################### 
####################### MODEL INCLUDING GENETIC DIVERSITY 
####################### 

#Test oversdispersion
fit4 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + (1|Block), data = data, family = "poisson")
blmeco::dispersion_glmer(fit4)
## [1] 5.69289
mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
blmeco::dispersion_glmer(mod)
## [1] 1.075746
#Convergence issue 
max(abs(with(mod@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001
## [1] 0.07273844
#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100), 
                       Generation_Eggs = rep(seq(1,6, length = 100),each=3),
                       Block = "Block4", 
                       Estimates = NA)
filldata$Estimates <- predict(fit4, newdata=filldata, type = "response")


#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates, 
                               colour = Genetic_Diversity), size = 1.3) +
  geom_point(data = data, 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs*Genetic_Diversity + 
#              gen_square*Genetic_Diversity, 
#            data = data[data$Generation_Eggs>=2,])
# fit2_test <- lm(log(Nb_adults+1)~ Generation_Eggs+Genetic_Diversity + 
#              gen_square, 
#            data = data[data$Generation_Eggs>=2,])

mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
mod1 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity + 
               (1|Block) + (1|obs), data = data, family = "poisson")
anova(mod,mod1, test = "Chisq")
## Data: data
## Models:
## mod1: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity + 
## mod1:     (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity + 
## mod:     (1 | Block) + (1 | obs)
##      npar    AIC  BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod1    9 5541.3 5579 -2761.7   5523.3                         
## mod    17 5527.0 5598 -2746.5   5493.0 30.365  8  0.0001822 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Signifcant interaction 
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##   5527.0   5598.0  -2746.5   5493.0      466 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.12632 -0.04460  0.00513  0.06208  0.29142 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.58420  0.7643  
##  Block  (Intercept) 0.07837  0.2799  
## Number of obs: 483, groups:  obs, 483; Block, 3
## 
## Fixed effects:
##                                                                       Estimate
## (Intercept)                                                          -1.987774
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         10.850818
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         -5.143118
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.941456
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         -0.058946
## Genetic_DiversityMedium                                              -0.598092
## Genetic_DiversityLow                                                 -0.917922
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  1.230996
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium -0.776495
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  0.157397
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium -0.010610
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     1.729870
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow    -0.976354
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     0.171344
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow    -0.009432
##                                                                      Std. Error
## (Intercept)                                                            1.224862
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           1.974563
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                           1.029881
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           0.213080
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                           0.015165
## Genetic_DiversityMedium                                                1.822946
## Genetic_DiversityLow                                                   1.616291
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   2.977179
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   1.559898
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.324231
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.023172
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      2.630789
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      1.374871
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.285225
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.020358
##                                                                      z value
## (Intercept)                                                           -1.623
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           5.495
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          -4.994
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           4.418
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          -3.887
## Genetic_DiversityMedium                                               -0.328
## Genetic_DiversityLow                                                  -0.568
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   0.413
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  -0.498
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.485
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  -0.458
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      0.658
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     -0.710
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.601
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     -0.463
##                                                                      Pr(>|z|)
## (Intercept)                                                          0.104620
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         3.90e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         5.92e-07
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         9.95e-06
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         0.000101
## Genetic_DiversityMedium                                              0.742843
## Genetic_DiversityLow                                                 0.570090
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.679257
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.618635
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.627359
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.647046
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow    0.510829
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow    0.477616
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow    0.548017
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow    0.643123
##                                                                         
## (Intercept)                                                             
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         ***
## Genetic_DiversityMedium                                                 
## Genetic_DiversityLow                                                    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 1.3357 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
summary(mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##   5541.3   5579.0  -2761.7   5523.3      474 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.09399 -0.05085  0.01382  0.07441  0.23786 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.62340  0.7896  
##  Block  (Intercept) 0.07624  0.2761  
## Number of obs: 483, groups:  obs, 483; Block, 3
## 
## Fixed effects:
##                                               Estimate Std. Error z value
## (Intercept)                                  -2.164709   0.781262  -2.771
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 11.963215   1.252559   9.551
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.794296   0.658426  -8.800
## stats::poly(Generation_Eggs, 4, raw = TRUE)3  1.063434   0.137051   7.759
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.066353   0.009797  -6.773
## Genetic_DiversityMedium                      -0.510577   0.095339  -5.355
## Genetic_DiversityLow                         -0.641010   0.086094  -7.445
##                                              Pr(>|z|)    
## (Intercept)                                   0.00559 ** 
## stats::poly(Generation_Eggs, 4, raw = TRUE)1  < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2  < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 8.53e-15 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 1.26e-11 ***
## Genetic_DiversityMedium                      8.54e-08 ***
## Genetic_DiversityLow                         9.66e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                    (Intr) s::(G_E,4,r=TRUE)1 s::(G_E,4,r=TRUE)2
## s::(G_E,4,r=TRUE)1 -0.964                                      
## s::(G_E,4,r=TRUE)2  0.940 -0.993                               
## s::(G_E,4,r=TRUE)3 -0.914  0.976             -0.995            
## s::(G_E,4,r=TRUE)4  0.888 -0.957              0.984            
## Gntc_DvrstM        -0.063  0.009             -0.010            
## Gntc_DvrstL        -0.066  0.005             -0.006            
##                    s::(G_E,4,r=TRUE)3 s::(G_E,4,r=TRUE)4 Gnt_DM
## s::(G_E,4,r=TRUE)1                                             
## s::(G_E,4,r=TRUE)2                                             
## s::(G_E,4,r=TRUE)3                                             
## s::(G_E,4,r=TRUE)4 -0.997                                      
## Gntc_DvrstM         0.011             -0.011                   
## Gntc_DvrstL         0.007             -0.007              0.491
## convergence code: 0
## Model failed to converge with max|grad| = 0.338298 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE  df asymp.LCL asymp.UCL
##  High                4.57 0.203 Inf      4.08      5.05
##  Medium              3.94 0.212 Inf      3.44      4.45
##  Low                 3.71 0.197 Inf      3.24      4.18
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.624 0.182 Inf 3.434   0.0017 
##  High - Low       0.854 0.163 Inf 5.229   <.0001 
##  Medium - Low     0.231 0.177 Inf 1.306   0.3916 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
####################### 
####################### MODEL INCLUDING ONLY RESCUED POPULATIONS
####################### 
#Same but only for population not extinct
fit_fin <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity,
            data = data[data$Extinction_finalstatus=="no",], family = "poisson")

filldata$Estimates_Withoutextpop <- predict(fit_fin, newdata=filldata, type = "response")

PLOT_Pop_size <- ggplot2::ggplot(data = data) + 
  geom_line(data = filldata, aes(x = Generation_Eggs, 
                                y = Estimates_Withoutextpop, 
                               colour = Genetic_Diversity), size = 1.5, linetype = "longdash") +
  geom_point(data = data[data$Extinction_finalstatus=="no",], 
                      aes(x = Generation_Eggs, 
                                y = Nb_adults, 
                                colour = Genetic_Diversity),
             position = position_dodge(0.4), size = 1.4, alpha = 0.5) +
  ylab("Population size") +
  labs(color="Genetic diversity") +
  xlab("Generation") +
  theme_LO
PLOT_Pop_size

mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + 
               (1|Block) + (1|obs), 
             data = data[data$Extinction_finalstatus=="no",], family = "poisson")
mod1 <- lme4::glmer(Nb_adults~poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity + 
               (1|Block) + (1|obs), 
              data = data[data$Extinction_finalstatus=="no",], family = "poisson")
anova(mod,mod1, test = "Chisq")
## Data: data[data$Extinction_finalstatus == "no", ]
## Models:
## mod1: Nb_adults ~ poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity + 
## mod1:     (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity + 
## mod:     (1 | Block) + (1 | obs)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## mod1    9 4759.9 4796.3 -2370.9   4741.9                       
## mod    17 4756.3 4825.1 -2361.1   4722.3 19.599  8    0.01196 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(mod)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +  
##     (1 | Block) + (1 | obs)
##    Data: data[data$Extinction_finalstatus == "no", ]
## 
##      AIC      BIC   logLik deviance df.resid 
##   4756.3   4825.1  -2361.1   4722.3      407 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.57063 -0.04879  0.00393  0.06967  0.29650 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  obs    (Intercept) 0.36624  0.6052  
##  Block  (Intercept) 0.02456  0.1567  
## Number of obs: 424, groups:  obs, 424; Block, 3
## 
## Fixed effects:
##                                                                        Estimate
## (Intercept)                                                          -1.9250178
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         10.7535583
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         -5.0922588
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.9316535
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         -0.0583188
## Genetic_DiversityMedium                                               1.1892488
## Genetic_DiversityLow                                                  0.2794991
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium -2.0341925
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  1.0112012
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.2039827
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  0.0137179
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow    -0.4034287
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     0.0777534
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow    -0.0101264
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     0.0006044
##                                                                      Std. Error
## (Intercept)                                                           0.9705855
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                          1.5722049
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          0.8204596
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                          0.1698215
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          0.0120893
## Genetic_DiversityMedium                                               1.5049889
## Genetic_DiversityLow                                                  1.3816159
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  2.4448787
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium  1.2751857
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  0.2640143
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium  0.0188078
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     2.2495506
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow     1.1751948
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     0.2435505
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow     0.0173582
##                                                                      z value
## (Intercept)                                                           -1.983
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                           6.840
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                          -6.207
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                           5.486
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                          -4.824
## Genetic_DiversityMedium                                                0.790
## Genetic_DiversityLow                                                   0.202
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium  -0.832
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   0.793
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium  -0.773
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.729
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow     -0.179
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      0.066
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow     -0.042
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.035
##                                                                      Pr(>|z|)
## (Intercept)                                                            0.0473
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         7.93e-12
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         5.41e-10
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         4.11e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         1.41e-06
## Genetic_DiversityMedium                                                0.4294
## Genetic_DiversityLow                                                   0.8397
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium   0.4054
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium   0.4278
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium   0.4397
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium   0.4658
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow      0.8577
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow      0.9472
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow      0.9668
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow      0.9722
##                                                                         
## (Intercept)                                                          *  
## stats::poly(Generation_Eggs, 4, raw = TRUE)1                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3                         ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4                         ***
## Genetic_DiversityMedium                                                 
## Genetic_DiversityLow                                                    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium    
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow       
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 3.25364 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?
emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE  df asymp.LCL asymp.UCL
##  High                4.53 0.133 Inf      4.21      4.85
##  Medium              4.30 0.151 Inf      3.94      4.66
##  Low                 4.01 0.137 Inf      3.68      4.33
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.230 0.153 Inf 1.507   0.2874 
##  High - Low       0.523 0.140 Inf 3.746   0.0005 
##  Medium - Low     0.293 0.157 Inf 1.861   0.1503 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
#Signifcant interaction 

3.4 Evolution of fitness over time

############ 
###### Distribution
############
hist(data_evolution_Lambda$Lambda, breaks=50)

# Square root
msqrt <- lme4::lmer(sqrt(Lambda) ~ Genetic_Diversity*Phenotyping  + 
                Block + (1|Population), 
              data = data_evolution_Lambda)
summary(msqrt)
## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(Lambda) ~ Genetic_Diversity * Phenotyping + Block + (1 |  
##     Population)
##    Data: data_evolution_Lambda
## 
## REML criterion at convergence: 252.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4979 -0.4737  0.0672  0.5620  2.4161 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Population (Intercept) 0.0295   0.1718  
##  Residual               0.1040   0.3225  
## Number of obs: 301, groups:  Population, 88
## 
## Fixed effects:
##                                           Estimate Std. Error t value
## (Intercept)                               1.779696   0.081683  21.788
## Genetic_DiversityMedium                  -0.171161   0.103522  -1.653
## Genetic_DiversityLow                     -0.150003   0.094035  -1.595
## PhenotypingFinal                         -0.243939   0.070714  -3.450
## BlockBlock4                               0.102764   0.063883   1.609
## BlockBlock5                               0.003821   0.070206   0.054
## Genetic_DiversityMedium:PhenotypingFinal -0.251050   0.110515  -2.272
## Genetic_DiversityLow:PhenotypingFinal    -0.154999   0.101256  -1.531
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM Gnt_DL PhntyF BlckB4 BlckB5 G_DM:P
## Gntc_DvrstM -0.600                                          
## Gntc_DvrstL -0.655  0.506                                   
## PhntypngFnl -0.685  0.534  0.588                            
## BlockBlock4 -0.462  0.054  0.041  0.028                     
## BlockBlock5 -0.424  0.006  0.009  0.002  0.488              
## Gntc_DvM:PF  0.428 -0.745 -0.377 -0.641 -0.028  0.050       
## Gntc_DvL:PF  0.444 -0.373 -0.733 -0.699  0.012  0.088  0.453
# log
mlog <- lme4::lmer(log(Lambda+1) ~ Genetic_Diversity*Phenotyping  +
               Block + (1|Population),
              data=data_evolution_Lambda)

# Normal
mnorm <- lme4::lmer(Lambda ~ Genetic_Diversity*Phenotyping  + 
                Block + (1|Population),
              data=data_evolution_Lambda)

## Compare AIC
AIC(mlog, msqrt,mnorm)
##       df      AIC
## mlog  10 207.0029
## msqrt 10 272.0671
## mnorm 10 848.9685
#Square root a better fits to the data 
#But we compare different values: Growth rate and Growth rate+1

#Distribution Residuals
hist(residuals(mlog),col="yellow",freq=F)

hist(residuals(msqrt),col="yellow",freq=F)  #Seems the better fit

hist(residuals(mnorm),col="yellow",freq=F)

#Test normality
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mlog), "norm"))
g1$kstest 
## 1-mle-norm 
## "rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(msqrt), "norm"))
g1$kstest 
##     1-mle-norm 
## "not rejected"
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest
##     1-mle-norm 
## "not rejected"
# QQ plot
# plot(mlog, which = 2) 
# plot(msqrt, which = 2) #Seems the better fit!
# plot(mnorm, which = 2) #Seems the better fit but very high AIC 
qqnorm(resid(mlog))
qqline(resid(mlog))

qqnorm(resid(msqrt))
qqline(resid(msqrt))

qqnorm(resid(mnorm))
qqline(resid(mnorm)) #best fit

############ 
###### Model selection
############

m0 <- lme4::lmer(Lambda ~ Genetic_Diversity*Phenotyping  + 
                Block + (1|Population) ,
              data=data_evolution_Lambda)


m1 <- lme4::lmer(Lambda ~ Genetic_Diversity+Phenotyping  + 
                Block + (1|Population) ,
              data=data_evolution_Lambda)

anova(m0,m1,test="Chisq") #Remove interaction
## Data: data_evolution_Lambda
## Models:
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
## m0: Lambda ~ Genetic_Diversity * Phenotyping + Block + (1 | Population)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## m1    8 832.93 862.59 -408.46   816.93                     
## m0   10 834.25 871.32 -407.12   814.25 2.6798  2     0.2619
#Same result with: 
# m1 <- lmer(Lambda ~ Genetic_Diversity+Phenotyping  + 
#                 (1|Block) + (1|Population) + (1|ID_Rep:Population),
#               data=data_evolution_Lambda)
# OR
# m1 <- lmer(Lambda ~ Genetic_Diversity+Phenotyping  + 
#                (1|Block) + (1|Population) ,
#               data=data_evolution_Lambda)

m2 <- lme4::lmer(Lambda ~ Phenotyping  + 
                Block + (1|Population),
              data=data_evolution_Lambda)
m3 <- lme4::lmer(Lambda ~ Genetic_Diversity  + 
                Block + (1|Population),
              data=data_evolution_Lambda)

anova(m1,m2,test="Chisq") #Keep phenotyping
## Data: data_evolution_Lambda
## Models:
## m2: Lambda ~ Phenotyping + Block + (1 | Population)
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m2    6 855.30 877.54 -421.65   843.30                         
## m1    8 832.93 862.59 -408.46   816.93 26.369  2   1.88e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(m1,m3,test="Chisq") #Keep genetic diversity
## Data: data_evolution_Lambda
## Models:
## m3: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## m3    7 896.49 922.44 -441.24   882.49                         
## m1    8 832.93 862.59 -408.46   816.93 65.558  1  5.642e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
mod_final <- lme4::lmer(Lambda ~ Genetic_Diversity+Phenotyping  + 
                Block + (1|Population),
              data=data_evolution_Lambda)

emmeans::emmeans(mod_final, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE   df lower.CL upper.CL
##  High                3.05 0.120 69.4     2.76     3.34
##  Medium              2.14 0.143 87.6     1.79     2.48
##  Low                 2.32 0.123 99.8     2.02     2.62
## 
## Results are averaged over the levels of: Phenotyping, Block 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE   df t.ratio p.value
##  High - Medium    0.916 0.184 72.8  4.988  <.0001 
##  High - Low       0.734 0.170 78.2  4.316  0.0001 
##  Medium - Low    -0.182 0.187 87.5 -0.972  0.5963 
## 
## Results are averaged over the levels of: Phenotyping, Block 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: tukey method for comparing a family of 3 estimates
emmeans::emmeans(mod_final, list(pairwise ~ Phenotyping), adjust = "tukey") #
## $`emmeans of Phenotyping`
##  Phenotyping emmean     SE  df lower.CL upper.CL
##  Initial       3.00 0.1072 254     2.76     3.24
##  Final         2.01 0.0837 109     1.82     2.20
## 
## Results are averaged over the levels of: Genetic_Diversity, Block 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 2 estimates 
## 
## $`pairwise differences of Phenotyping`
##  contrast        estimate    SE  df t.ratio p.value
##  Initial - Final    0.993 0.117 250 8.474   <.0001 
## 
## Results are averaged over the levels of: Genetic_Diversity, Block 
## Degrees-of-freedom method: kenward-roger
#

3.5 Fitness G6

m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    Data: data_phenotyping
## 
##      AIC      BIC   logLik deviance df.resid 
##   2154.4   2171.3  -1072.2   2144.4      212 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1441 -0.1230  0.0254  0.1360  0.3923 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  ID_Rep:Population (Intercept) 0.2041   0.4518  
##  Population        (Intercept) 0.2807   0.5298  
## Number of obs: 217, groups:  ID_Rep:Population, 217; Population, 77
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               4.2792     0.1105  38.715  < 2e-16 ***
## Genetic_DiversityMedium  -0.7386     0.1812  -4.076 4.59e-05 ***
## Genetic_DiversityLow     -0.5250     0.1664  -3.156   0.0016 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM
## Gntc_DvrstM -0.611       
## Gntc_DvrstL -0.667  0.413
# Equivalent to:
# m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population/ID_Rep),
#                 family = "poisson", data = data_5week)


#dispersion_lme4::glmer(m0) #dispersion ok 
#Note: (1|School) + (1|Class:School) is equivalent to (1|School/Class)

m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Population) + (1|ID_Rep:Population),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") # difference
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Population) + (1 | ID_Rep:Population)
## m0: Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)    
## m1    3 2167.5 2177.6 -1080.7   2161.5                        
## m0    5 2154.4 2171.3 -1072.2   2144.4 17.04  2  0.0001994 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(m0)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
##    Data: data_phenotyping
## 
##      AIC      BIC   logLik deviance df.resid 
##   2154.4   2171.3  -1072.2   2144.4      212 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1441 -0.1230  0.0254  0.1360  0.3923 
## 
## Random effects:
##  Groups            Name        Variance Std.Dev.
##  ID_Rep:Population (Intercept) 0.2041   0.4518  
##  Population        (Intercept) 0.2807   0.5298  
## Number of obs: 217, groups:  ID_Rep:Population, 217; Population, 77
## 
## Fixed effects:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               4.2792     0.1105  38.715  < 2e-16 ***
## Genetic_DiversityMedium  -0.7386     0.1812  -4.076 4.59e-05 ***
## Genetic_DiversityLow     -0.5250     0.1664  -3.156   0.0016 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Gnt_DM
## Gntc_DvrstM -0.611       
## Gntc_DvrstL -0.667  0.413
############ 
###### Posthoc
############ 
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean    SE  df asymp.LCL asymp.UCL
##  High                4.28 0.111 Inf      4.02      4.54
##  Medium              3.54 0.143 Inf      3.20      3.88
##  Low                 3.75 0.124 Inf      3.46      4.05
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate    SE  df z.ratio p.value
##  High - Medium    0.739 0.181 Inf  4.076  0.0001 
##  High - Low       0.525 0.166 Inf  3.156  0.0046 
##  Medium - Low    -0.214 0.189 Inf -1.132  0.4943 
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates
#test double nested
m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Genetic_Diversity/Population/ID_Rep),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Genetic_Diversity/Population/ID_Rep),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") 
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/Population/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/Population/ID_Rep)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## m1    4 2160.9 2174.4 -1076.4   2152.9                       
## m0    6 2156.4 2176.7 -1072.2   2144.4 8.4441  2    0.01467 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
m0 <- lme4::glmer(Adults ~ Genetic_Diversity  + (1|Genetic_Diversity/ID_Rep),
                family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1  + (1|Genetic_Diversity/ID_Rep),
                family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") 
## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/ID_Rep)
##    npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## m1    3 2189.9 2200.0 -1091.9   2183.9                        
## m0    5 2183.0 2199.9 -1086.5   2173.0 10.858  2   0.004387 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

3.6 Fitness G2

m0 <- glm(Nb_adults ~ Genetic_Diversity + Block,  family = "poisson", data = data_G2)
summary(m0)
## 
## Call:
## glm(formula = Nb_adults ~ Genetic_Diversity + Block, family = "poisson", 
##     data = data_G2)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -18.7057   -2.7621    0.2126    3.3983   11.4454  
## 
## Coefficients:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              5.73890    0.01497 383.312  < 2e-16 ***
## Genetic_DiversityMedium -0.19408    0.01628 -11.920  < 2e-16 ***
## Genetic_DiversityLow    -0.19278    0.01460 -13.205  < 2e-16 ***
## BlockBlock4              0.10767    0.01579   6.817  9.3e-12 ***
## BlockBlock5              0.17743    0.01600  11.088  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 2379.6  on 83  degrees of freedom
## Residual deviance: 2027.9  on 79  degrees of freedom
## AIC: 2667.2
## 
## Number of Fisher Scoring iterations: 4
m1<- glm(Nb_adults ~ Block,  family = "poisson", data = data_G2)
anova(m0,m1,test="Chisq") # difference
## Analysis of Deviance Table
## 
## Model 1: Nb_adults ~ Genetic_Diversity + Block
## Model 2: Nb_adults ~ Block
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1        79     2027.9                          
## 2        81     2241.4 -2  -213.52 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean      SE  df asymp.LCL asymp.UCL
##  High               5.834 0.01051 Inf     5.809     5.859
##  Medium             5.640 0.01243 Inf     5.610     5.670
##  Low                5.641 0.01022 Inf     5.617     5.666
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate     SE  df z.ratio p.value
##  High - Medium   0.1941 0.0163 Inf 11.920  <.0001 
##  High - Low      0.1928 0.0146 Inf 13.205  <.0001 
##  Medium - Low   -0.0013 0.0161 Inf -0.081  0.9964 
## 
## Results are averaged over the levels of: Block 
## Results are given on the log (not the response) scale. 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.7 Remaining He difference

Rmisc::summarySE(data_G2,
          measurevar = c("He_remain"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N He_remain          sd          se          ci
## 1              High 27 0.9986008 0.000000000 0.000000000 0.000000000
## 2            Medium 23 0.6791246 0.006059307 0.001263453 0.002620241
## 3               Low 34 0.5441887 0.012687849 0.002175948 0.004427000
mfull <- lm(He_remain ~ Genetic_Diversity, data=data_G2)
mless <- lm(He_remain ~ 1, data=data_G2)
anova(mfull,mless) #Keep effect
## Analysis of Variance Table
## 
## Model 1: He_remain ~ Genetic_Diversity
## Model 2: He_remain ~ 1
##   Res.Df    RSS Df Sum of Sq     F    Pr(>F)    
## 1     81 0.0061                                 
## 2     83 3.1868 -2   -3.1807 21048 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean       SE df lower.CL upper.CL
##  High              0.9986 0.001673 81   0.9945   1.0027
##  Medium            0.6791 0.001812 81   0.6747   0.6835
##  Low               0.5442 0.001491 81   0.5406   0.5478
## 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate      SE df t.ratio p.value
##  High - Medium    0.319 0.00247 81 129.527 <.0001 
##  High - Low       0.454 0.00224 81 202.800 <.0001 
##  Medium - Low     0.135 0.00235 81  57.498 <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
Rmisc::summarySE(data_G2,
          measurevar = c("He_remain"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N He_remain          sd          se          ci
## 1              High 27 0.9986008 0.000000000 0.000000000 0.000000000
## 2            Medium 23 0.6791246 0.006059307 0.001263453 0.002620241
## 3               Low 34 0.5441887 0.012687849 0.002175948 0.004427000
#Calcul for end:
Rmisc::summarySE(data_G2[data_G2$Extinction_finalstatus!="yes",],
          measurevar = c("He_remain_end"),
          groupvars = c("Genetic_Diversity"),
          na.rm=TRUE)
##   Genetic_Diversity  N He_remain_end         sd          se          ci
## 1              High 27     0.9697741 0.01868272 0.003595491 0.007390637
## 2            Medium 18     0.6482803 0.02453741 0.005783522 0.012202164
## 3               Low 26     0.5146307 0.02759467 0.005411760 0.011145728
mfull <- lm(He_remain_end ~ Genetic_Diversity, data=data_G2[data_G2$Extinction_finalstatus!="yes",])
mless <- lm(He_remain_end ~ 1, data=data_G2[data_G2$Extinction_finalstatus!="yes",])
anova(mfull,mless) #Keep effect
## Analysis of Variance Table
## 
## Model 1: He_remain_end ~ Genetic_Diversity
## Model 2: He_remain_end ~ 1
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1     68 0.03835                                  
## 2     70 2.91180 -2   -2.8735 2547.7 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
############ 
###### Posthoc
############ 
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") 
## $`emmeans of Genetic_Diversity`
##  Genetic_Diversity emmean      SE df lower.CL upper.CL
##  High               0.970 0.00457 68    0.959    0.981
##  Medium             0.648 0.00560 68    0.635    0.662
##  Low                0.515 0.00466 68    0.503    0.526
## 
## Confidence level used: 0.95 
## Conf-level adjustment: sidak method for 3 estimates 
## 
## $`pairwise differences of Genetic_Diversity`
##  contrast      estimate      SE df t.ratio p.value
##  High - Medium    0.321 0.00723 68 44.491  <.0001 
##  High - Low       0.455 0.00653 68 69.754  <.0001 
##  Medium - Low     0.134 0.00728 68 18.355  <.0001 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

3.8 Growth rate and He

# 
# #fit first degree polynomial equation:
# fit  <- lm(Lambda ~ HE, data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# #same fit as: fit5 <- lm(Lambda~HE+HE^2, data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# #same fit as: fit6 <- lm(Lambda~HE^2, data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# 
# #second degree
# fit2 <- lm(Lambda~poly(HE,2,raw=TRUE), data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# #third degree
# fit3 <- lm(Lambda~poly(HE,3,raw=TRUE), data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# #fourth degree
# fit4 <- lm(Lambda~poly(HE,4,raw=TRUE), data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# #Other fit 
# fit7 <- lm(Lambda~exp(HE), data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# fit8 <- lm(Lambda~log(HE), data = SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",])
# 
# #generate range of 50 numbers starting from 30 and ending at 160
# xx <- seq(0.3,1, length=100)
# plot(SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",]$HE,
#      SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",]$Lambda,pch=19,ylim=c(0,5))
# lines(xx, predict(fit, data.frame(HE=xx)), col="red")
# lines(xx, predict(fit2, data.frame(HE=xx)), col="green")
# lines(xx, predict(fit3, data.frame(HE=xx)), col="blue")
# lines(xx, predict(fit4, data.frame(HE=xx)), col="purple")
# # lines(xx, predict(fit5, data.frame(HE=xx)), col="black")
# # lines(xx, predict(fit6, data.frame(HE=xx)), col="yellow")
# lines(xx, predict(fit7, data.frame(HE=xx)), col="brown")
# lines(xx, predict(fit8, data.frame(HE=xx)), col="orange")
#